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, Abstract. The relativistic time evolution of multi-layer spherically symmetric shell 

' systems — consisting of infinitely thin shells separated by vacuum regions — is examined. 

Whenever two shells collide the evolution is continued with the assumption that 
the collision is totally transparent. The time evolution of various multi- layer shell 
O 1 systems — comprised by large number of shells thereby mimicking the behavior of 

a thick shell making it possible to study the formation of acoustic singularities — is 
analyzed numerically and compared in certain cases to the corresponding Newtonian 
time evolution. The analytic setup is chosen such that the developed code is capable of 
following the evolution even inside the black hole region. This, in particular, allowed 
us to investigate the mass inflation phenomenon in the chosen framework. 
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1. Introduction 

o 

y- —I ' 

Relativistic infinitely thin spherical shells play an important role in various dynamical 
contexts ranging from microscopic to astrophysical systems. For instance, by applying 
a charged shell as an electron model one may avoid the appearance of negative 
gravitational mass caused by the concentration of charge at the center [TJ [2j [3]. 
Using families of spherically symmetric thin shells instead of spherically symmetric 
continuous matter distributions reduces significantly the complexity of evolutionary 
problems as the dynamics of thin shells may be investigated by using various analogies 
from the description of the motion of a particle in a one- dimensional effective potential. 
Accordingly, the quantization of systems comprising thin shells is tractable [U |2j [5j [6] . 
Macroscopically stable quark-gluon matter can also be studied with a toy model in 
which relativistic shells and the MIT bag model are combined [7J. Collapsing dust 
shells can be used to probe stability or to study energetics of compact objects such 
as back holes or star models mimicking the properties of black holes [8] E2J [9] . Shells 
can be used to model matter ejection at certain phases of supernova explosions [10] or 
in modeling supernova remnants [TT]. More realistic radiating shell models can also 
be constructed [12] and with the help of these models, even the critical collapse may 
be investigated analytically |13j . With the help of infinitesimally thin shells one can 
construct exact solutions by gluing together spherically symmetric spacetime domains. 
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This way exotic models such as gravastars [HI [33] or wormholes [TBI [TBI [Till l2TJ l l2T l l2"iZ l l2"3"] 
may also be studied. Simple models of large-scale voids in galaxy distributions can also 
be constructed with the help of shells [241 1251 126] . The dynamics of spherical shells come 
into play in some cosmological models, such as higher dimensional brane cosmologies, 
in which it is assumed that our four-dimensional universe is merely a surface living in 
a higher dimensional spacetime [27]. Shells play a central role in the bubble inflation 
model of the early universe [281 1211 [301 [311 [32] . It is widely thought that by studying 
dynamics of shells, important phenomena such as the focusing singularity at the center 
[33] or the so-called acoustic singularity — which will be discussed later in section 4.1.— 
can also be studied. 

The basic equations governing the dynamics of thin shells can be derived in 
various ways, for instance by making use of the junction conditions of the metrics 
[SI EHl EEl EHl EHl SQl SH S21 S3] , by applying distributions [H SSI HE] or the variational 
Jacobi-Hamiltonian approach [ITJ SH HH [501 [51], or by deriving them by taking a limit 
of the evolution equations of the thick shells [52, 53J. It is important to mention that 
all of these methods produce the same set of basic equations. 

In describing the evolution of families of infinitesimally thin shells, the study of 
their crossing is essential. Nevertheless, much less has been done in this respect. For 
instance, the basic equations describing shell crossing have only been derived for some 
specific cases, such as totally transparent [HH EH [55] or totally inelastic collisions [561 E] ■ 
Nevertheless, most of the authors do not go beyond deriving the equations of motion 
for dust shells, or studying only the simplest possible analytic cases. This, in particular, 
means that almost no results are available for multi-layer shell systems with the generic 
equation of state (EOS). There are, apparently, only two papers in the literature which 
carry out a systematic study of the dynamics of more than two shells or consider the 
evolution of shell systems such that shell crossings are allowed. In [51] the dynamics 
of star clusters is studied, although considerations therein remain on the theoretical 
side, providing only the generic equations of motion in Schwarzschild time coordinates 
which, in particular, does not allow the study of the motion through the event horizon. 
In [55] the dynamics of two dust shells was investigated with the use of Kruskal-Szekeres 
coordinates — thereby the authors were able to follow the motion of the system below 
the horizon — and shell crossings were assumed to be totally transparent. 

Our main aim in this paper is to present some new results concerning the dynamics 
of multi-layer shell systems with the inclusion of of shell crossings. The corresponding 
dynamical investigations were carried out by using a C++ code [57] which made the 
study of the evolution of systems composed of a large number of shells and with a 
generic EOS possible. In section 2, some of the basics related to the analytic description 
of the motion of a single shell are recalled using Schwarzschild and ingoing Eddington- 
Finkelstein 'time' coordinates. For comparison, the corresponding Newtonian case is 
also discussed. In section 3, for the sake of simplicity, first only the interaction of two 
shells is considered, providing the balance equations for the case of totally transparent 
shell crossings. In section 4, the time evolution of systems comprising a large number of 
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concentric shells — in particular, the formation of acoustic singularities — is investigated. 
One of the main advantage of the method applied in this paper is related to the use of 
ingoing Eddington-Finkelstein coordinates, which allows us to evolve even very complex 
composite systems throughout the entire spacetime, including the black hole region. One 
of the most important outcomes of the investigations is that the phenomena of mass 
inflation could also be studied by making use of colliding thin shells. The paper is 
concluded by our final remarks. 

Throughout this paper, the geometrized units are used, with G = c = 1, and the 
abstract index notation of [58] will be applied with the additional use of uppercase Latin 
indices signifying quantities living on three-dimensional hyper surf aces. 



2. The dynamics of a single shell 

In this section the basic equations relevant for a spherically symmetric infinitely thin 
shell, bounded by two Schwarzschild vacuum regions, will be recalled. 



2.1. Equation of motion 

Consider now a single shell and assume that the metric of the Schwarzschild regions on 
the sides of the shell is given as 

d4 = -f±(r) dt 2 ± + / ± (r)" 1 dr 2 + r 2 (dtf 2 + sin^d^ 2 ), (2.1) 

were the indices ± signify the outer and inner regions, respectively, r stands for the 
area radius, which is assumed to be continuous across the shell, t_ and t + denote the 
Schwarzschild time coordinates in the corresponding spacetime regions, *& and ip are the 
standard spherical coordinates, while 

Mr) = 1 - ^ (2.2) 

r 

with mass parameters M±, respectively^. 

Denote by u a the four-velocity tangent to the timelike generators of the shell and 
by r the proper time along these timelike generators. The components of u a can be 
given as 

<-(§■£■*»)■ M> 

The induced metrics, h~ AB and h\ Bl on the mutual boundary of the two spacetime 
regions are assumed to coincide, and the metric on the shell Hab = h\ B = h AB , in the 
(r, cp) coordinates, reads 

h AB = diag(-l, r 2 (T),r 2 (r) sin 2 #) . (2.4) 

1 Note that the form of the metric (|2.1j) with the slightly more generic metric function f±(r) = 
1 — 2M±(r)/r is suitable to represent besides the Schwarzschild metric, that of the de Sitter spacetime 
in the vacuum case, or whenever electrovacuum spacetimes are also included the Reisner-Nordstrom 
de Sitter solutions also fit this form. 
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(*-.»•-) (U,r + ) 

Figure 1. The convention applied for fixing the sign of e n± is illustrated. The regions 
to be matched are indicated by the grey domains while £± denotes their 'isometric' 
boundaries in the preliminary spacetimes. 



As the four-velocity u a is of unit norm its components are not independent and, in virtue 
of (12. ip and f l2.3|) . the relation u a u a = — 1 implies that the relation 

/ ± ^V=/ ± +C?V (2-5) 



dr / \dr J 

holds from which, with the inclusion of the sign factor e^ ± = sign(/± • dt±/dr 
dt ± _ e i±v // ± (r) + (dr/dr)2 



(2.6) 



dr f ± (r) 

can be deduced. We shall return to the fixing of this sign factor later in subsection l2.5l 
Nevertheless, let us mention here only that the value of e*_ or e t+ is +1 in regions where 
t- or t + is a timelike coordinate, i.e. above the respective horizons. 

The unit normal n a to the shell, satisfying the orthogonality requirement u a n a = 0, 
can be given as 

n a ± = e n± (--, ^±,0,0j , (2.7) 

where the value of the sign factor e n± is chosen such that n a points outward as it is 
illustrated in figured! In the particular Schwarzschild case this implies that e n± = +1.§ 
Now, with the help of the induced metric, Hab-, an d the unite normal n a ± the 
extrinsic curvature tensors, at the boundaries £_ and £+, are given as 

K-ab = -jhA E h B F £ n± hEF ■ (2.8) 

2 Note that by allowing e n± to take the value ±1 within the very same setup, wormholes [HI [19] and 
other type of exotic spacetimes can be studied. However, in this paper attention will be restrected to 
conventional shells. 
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As already mentioned the equation of motion for the shell separating the 
Schwarzschild regions can be derived in various ways. One of the most frequently applied 
methods is based on Israel's thin-shell formalism pi] in which junction conditions 
are specified at the location of the shell. In particular, besides the relation Kab = 
h\ B = h~ AB , the discontinuity of the pertinent extrinsic curvature tensors, K\ B and 
K AB , is assumed to be related to the surface energy- momentum tensor of the shell 
Sab = diag(cx, V-r 2 , V-r 2 sin 2 6) — where a and V denote the surface energy density 
and the two-dimensional tangential pressure of the shell, respectively — via the Lanczos 
equation 

K\b ~ K AB = -8tt jsUfl - \h AB (h CD S CD ) } • (2.9) 

By making use of (12.71) and (I2.8P we get = e n ± r /±(dt±/dr). Then, in virtue of 
( 12. 9p . the component of this junction condition can be seen to take the form 




r)+(£) = 4 ^ ( 2 - 10 ) 

where r now signifies the radius at the location of the shell, while the sign factors s_ 
and s + are nothing but the products of e t± and e n± , i.e. s± = e t± e n± . Since e n± = +1, 
s± = e t± hereafter. 

Note that the importance of the appropriate treatment of the sign factors s_ and s + 
was emphasized in [36[ 07] . Regardless of which of the alternative methods (mentioned 
in the introduction) is applied, the reasonings always ends up with the equation of 
motion 

\dr J \ m r / r V 2r / 

where m c = M_ denotes the Schwarzschild mass parameter of the central region, 
m T = Aitar 2 is the 'rest mass' of the shell representing the surface internal energy 
associated with the tangential motion of particles, while m g = M + — M_ is the 
'gravitational mass' of the shell. Accordingly, m c and m g are constants but in general 
m r is a function of the radius, and it is constant only in the particular case of dust shells 
with zero pressure. 

To determine the r-dependence of m r we need an additional equation which can 
be derived by making use of the conservation law D a Sab = 0, where Da denotes 
the covariant derivative associated with the metric Hab- 111 this particular case the 
conservation law takes the form [39 rl 



^K)=^«. (2-12) 

3 Equation ([2.12[) may also be derived by making use of the two algebraically independent components 
of (1231). 
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By introducing the area radius r, instead of r, as our independent variable ( I2.12p can 
be put into the form 

r ^- = - 2 (a + V). (2.13) 
dr 

This is the point where the EOS of the shell comes into play. With the help of an 
EOS of the form V = V(a) the functional form of a = a(r) may be determined, 
which, in turn, gives us the functional form of m r (r), as well. Since (I2.13P is a separable 
differential equation an implicit solution to it — provided that the EOS is regular enough 
to guarantee the above integral to exist — can be written as 

(2.14) 

with integration constant cr = cr(r ). 

Once m T (r) is known from (12.111) . along with suitably chosen initial conditions, can 
be used to determine the motion as a function of the proper time r. 

2.2. Equation of state 

Although the developed C++ code allows the use of basically any kind of EOS in the 
numerical investigations covered by this paper, only the homogeneous linear EOS is 
applied. It is also important to keep in mind that in specifying an EOS some additional 
restrictions always need to be taken into account. For instance, the use of a suitable 
energy condition is essential. The most appropriate one is the so-called dominant energy 
condition (DEC) guaranteeing that solutions to dynamical problems respect the concept 
of causality. For a selected type of infinitesimally thin shell DEC can be seen to hold 
whenever \V\ < a, which, in particular, means that a is non-negative (see, e.g., [91 159]). 
To avoid dynamical instabilities we shall also assume that the square of the speed 
of sound, Cg = dV/da, is non-negative, and, to be compatible with the conventional 
concept of relativity, that cl is less than or equal to the square of the speed of light. 

As mentioned above, for the shake of simplicity, in this paper considerations will 
be restricted to the simplest possible functional form, i.e. to a homogeneous linear EOS 
V(cr) = wa in which case (12.141) takes the form 

t \ -2(l+to) 

a(r) = a (-) . (2.15) 

DEC, along with hydrodynamical stability on the surface, can be seen to hold whenever 
w is chosen from the interval < w < 1. Note also that w = corresponds to the dust 
case. 
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2.3. Initial conditions 

In applying Israel's junction conditions one may start by fixing the geometrical 
properties of the spacetimes on the respective sides of the shell, and then solving the 
junction conditions for suitable rest mass and initial velocity values making the matching 
possible. This could be referred to as a 'mathematician approach', while in a 'physicist 
approach' one would start by fixing the rest mass and initial velocity of the shell and 
then try to determine suitable geometrical parameters of the respective side spacetime 
regions. The latter approach is applied below. 

Accordingly, in describing the motion of a shell we regard m c as the environmental 
parameter and mo, ro, vq as initial parameters of the shell at To. Here mo = m r (ro), ro = 
r(r ), Vq = r(r ), and the over-dot denotes derivative with respect to the proper time 
t. As mentioned above, to be compatible with DEC we shall assume that mo > 0. In 
virtue of equation (12. lip , or (I2.10p . whenever the junction is possible the gravitational 
mass, m g , depends only on the kinetic energy and it may be expressed in terms of the 
initial data as 



m p 



m 



mo 
2r 



1 - + vl 



r 



(2.16) 



In order to avoid the gravitational mass to becoming complex the inequality 
2m 



v 2 Q > 



r 



- 1 



(2.17) 



has to hold. 

In virtue of ( 12. 16ft m g may be negative. Note also that since a is required to be 
positive e f _ = — 1 and e i+ = +1 cannot occur. As we shall see below, (I2.25|) excludes the 
possibility of having both e t _ and e t+ be negative. Thus, in the remaining two cases, the 
sufficient conditions ensuring the existence of the desired type of matching are given as 







sufficient conditions 


gravitational mass 


+ 


+ 


v 2 >A 


m g > 


+ 




v 2 <A 


m g > 0, if Vq > B\ m g < 0, elsewhere 



where 



.4 



B 



VTLr 



+ 



4r 2 



+ 



2m c 
2m c 



- 1 . 



(2.18) 
(2.19) 



2-4- Characterizing the dynamics 

For dust shells the values of m g and m r — the latter is also constant for dust shells — 
characterize the motion in the following way. The system is said to be gravitationally 
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bound whenever there exists a finite radius, r max , such that the velocity vanishes at 
r max . In general, the value of r max coincides with the positive root of the right-hand side 
of ( 12.1 ip which, whenever m g < m r reads as@ 



1_ ^J Uc + ^ + ^ + mcmg + y ) ■ (2.20) 

If m g = m r the kinetic energy becomes zero exactly at the spatial infinity, i.e. 
whenever r max = oo, and the associated motion is usually referred to as 'marginally 
bound'. For m g > m r the motion of the shell is not restricted, and, in virtue of ( I2.16p . 
at the spatial infinity, i.e. in the tq — > oo limit, the relation m g = m r a/1 + holds. 

Although the solution to ( 12. lip is, in general, complicated, for dust shells it can be 
given in closed form (see, e.g, [S2]). For example, in the special case of marginally bound 
motion, with m„ = m r , the pertinent solution can be given by the implicit relation 



(4rm c + 2rm r — m!?)\/8rm c + 4rm r + ml 
r(r) = 



2Am^ml 



(2.21) 



where r denotes the initial location of the shell at r = 0. Setting r = the value of 
— r(r) becomes equal to the proper time duration meanwhile the shell collapses from 
radius r to the symmetry center. Specializing even further, for a 'self-gravitating' shell 
with a Minkowski interior, i.e. with m c = 0, the proper time that is needed for such a 
shell of radius r to undergo a complete gravitational collapse can be given as 



m r y/4r + m r ( r ^/m r ] 

r c {r) = — + ~ — ■ 2.22 

6 3 \ ^nrir 2 



It is worth keeping in mind that the above analytic expressions were derived by making 
use of the assumption that the motion is marginally bound which, in particular, means 
that the velocity of the shell at r should be as if the shell started to move towards the 
center from rest at spatial infinity. 

Figure 2 shows some numeric examples for dust shells. For shells with non- 
zero pressure, the situation is more complex because pressure may, and in fact does, 
significantly alter the motion of the particles. By investigating the case of a homogeneous 
linear EOS it was justified, in contrast to the dust case, that such a shell may be in 
equilibrium, although the pertinent equilibrium was found to be unstable [60, 40J. A 
more complicated EOS such as polytrop can only be studied numerically, and — as was 
justified by our numerical experiences — it is possible to construct shells which oscillate 
in a bounded region or which are in stable equilibrium. 



2.5. The Schwarzschild time and the Eddington-Finkelstein null coordinates 

Up to this point all the derivatives in the equations relevant to the evolution of the 
investigated shell were expressed with respect to the proper time associated with the 

4 It follows from ()2.20|) that r max > Rs — 2(m c + m g ) as mentioned before. 
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shell. In many practical cases it turns out to be necessary — as it will clearly be 
demonstrated in the following sections concerning the collision of shells — to know what 
is the functional relation of this proper time, e.g., to the Schwarzschild time coordinates 
defined on the inner and outer sides of the shell. We would like to emphasize that in 
the applied formalism only the area radius, r, is required to be a continuous — although, 
not necessarily monotonic — function through the shells. Accordingly, in general the 
Schwarzschild time coordinates defined on the sides need not match continuously. 
Therefore, we have to determine the functional relation of the proper time to both 
the inner and outer Schwarzschild time coordinates separately. 

In doing so recall that the derivative (di±/dr) has been given by ( 12.61) . Note, 
however, that the value of e t± — for its definition see subsection l2.U — is undetermined yet, 
although it is uniquely determined in regions where the Schwarzschild time coordinates 
increase in the future direction, i.e. whenever the radius is larger than both of the 
Schwarzschild-radii, where dt±/dr > and f± > 0, the use of e i± = +1 in (12.61) is 
required. 

Nevertheless, as it was shown in [39] the derivative (dt±/dr) can always be 
determined uniquely. More concretely, by making use of the r — r component of (12.91) 
it can be shown — for details see, e.g., the part of section 2 in |39j between equations 
(2.21) and (2.34) — that regardless of the location of the shell, whenever it moves in a 
Schwarzschild spacetime 

dt_ / 2m c \ 1 / m„ m r \ , 

17 = V-^) U + 2f)' (2 ' 23) 

dt + ( ^ 2(m c + m g )\ 1 f m s m 



1 _ zx^L^JkL _ _i 1 (2.24) 

dr \ r \m r 2r 



or equivalently, the relations £t±yj±(^) + (dr/dr) 2 = m g /m r =p m r /2r are always 
satisfied everywhere. This, in particular, implies that 

e t± = sign(m g /m r =)= m r /2r) . (2.25) 

It follows from (12.251) that either e t _ or e t+ changes sign depending on the sign of 
m g . This sign change occurs at the 'critical radius' r = m 2 /(2|m g |).§ 

Once the initial values for t_ and t + are specified, these equations determine the 
desired relations between the proper time and the Schwarzschild time coordinates in the 
inner and outer spacetime regions, respectively. The differences of the right-hand sides 
of (I2.23P and (12.241) make it clear that, apart from very exceptional cases, even though 
the initial values for t- and t + are chosen to coincide, they will necessarily differ latter, 
i.e. they need not match continuously as was indicated above. 

The Schwarzschild metric has a coordinate singularity at the event horizon, i.e. 
at r = 2m c in the inner region and at r = 2(m c + m g ) in the outer region. The 

5 Note that whenever V > the value of r is unique due to the monotonicity of the functions 
m r = m r (r). It can also be checked that unless this change of the sign occurs, the worldsheets 
representing the evolving shells cannot be of class C 2 at r = r. 
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m c =0, m g =m r =1 , v(~)=0 m c =1 , m g =m r =1 , v(°o)=0 




time time 



Figure 2. In the left panel the motion of a single dust shell, starting at r = 10, in 
the particular case with m c = 0, m t = m g = 1 is shown with respect to the proper 
time, the interior Minkowski time and the exterior Schwarzschild time. In the last 
case, as is indicated by the pertinent plot, the shell can only reach the event horizon 
asymptotically, while the time of complete collapse is r c = (19\/41 + l)/6 ~ 20.44 
or f c = (llyil — l)/3 ~ 23.14 with respect to the proper or inner Minkowski time, 
respectively. In the right panel the motion of a similar shell with m c = m r = m g = 1 is 
indicated. In this case the inner and outer horizons denote the horizons at the central 
and exterior Schwarzschild regions, respectively. 



Schwarzschild time coordinate becomes infinite while approaching the horizon (see for 
example figure 2); thereby, it does not allow the proper description of the motion through 
the horizon. It is possible to overcome this technical difficulty by applying ingoing 
Eddington-Finkelstein coordinates covering both the domain of outer communication 
and the black hole regions simultaneously. The ingoing Eddington-Finkelstein null 
coordinate is defined as v± = t± + r*, where r* denotes the 'tortoise coordinate' 
determined by the relation dr*/dr = f± (r). Accordingly, 

dv± _ dt± | dr* dr _ e t±y / ' f±(r) + (dr/dr) 2 + dr/dr 
dr dr dr dr f±( r ) 

where e t± is given by (12.25p . 

Note that in the particular case of a dust shell with m c = and m g = m r , by 
making use of (12.231) and (I2.22p . the time needed for a complete gravitational collapse 
may be given in terms of the interior Minkowski time as 



\/4r + m r ( , r \ m r , _ s 

*c(r ) = ^—i -(V^+^) 2.27 



In concluding this section we would like to emphasize that in describing the motion 
of a single shell — determined by ( 12. lip — only the mass parameters of the interior and 
exterior Schwarzschild spacetime regions, along with the mass of the shell (determined 
by the surface energy density and the pertinent EOS of the shell), are relevant. In 
particular, as long as there is no collision between a selected and the surrounding shells, 
the motion of the selected one depends only on the mass parameters of the interior 
and exterior Schwarzschild regions Note also that, likewise in the Newtonian case, the 
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mass of the interior shells (if they exist at all) comes into play only via the value of the 
central mass m c . However, as the gravitational mass m g may not be positive, the general 
relativistic motion of the shell may differ significantly from that of a shell, possessing 
the same m c and m r parameters, in the Newtonian theory. 



2. 6. The motion of shells in the Newtonian case 

In this section, for the sake of completeness and for comparison, a brief review of the 
Newtonian framework will be given first. Then, the equations of motion will be derived 
for fluid shells. Throughout this subsection the basic quantities in the Newtonian 
framework will be signified by uppercase letters corresponding to the ones introduced 
in the previous sections for fully relativistic systems. 

Note first that in the Newtonian theory the inertial mass and the gravitational 
mass of a shell are not distinguished. This common mass of the shell — which, in the 
Newtonian description, is independent of time and the EOS — will be denoted by M S h e ii- 
The other considerable simplification characterizing the Newtonian theory comes from 
the use of absolute time, denoted by T. 

Likewise in the fully relativistic case the equation of motion of the shells in the 
Newtonian framework can be derived in various ways. In the case of dust shells, i.e. for 
shells with zero pressure, a good review can be found in |61j . According to the argument 
outlined therein the basic equation is nothing but a balance equation of force per unit 
mass and for the dust case it possesses the form of (12.281) below with P = 0. 

In generalizing this result to the fluid shells, i.e. in determining the functional 
form of the last term on the right-hand side of (12.281) with ? ^ 0, we need to find the 
appropriate force term representing the contribution of the pressure to the Newtonian 
balance equation. In identifying it let us consider an elementary 'square shaped' surface 
element of the shell with sides R A0, where A<p denotes the viewing angle of the sides 
from the center. The mass of this elementary piece is AM s h e ii = £-R 2 (A0) 2 , where 
£ = Af s hcii/(47ri? 2 ) is the surface mass density. There are four elementary forces exerted 
at the sides of this square-shaped surface element. The size of these elementary forces is 
AF = P RA(j), where P denotes the 'two-dimensional' pressure of the shell. Since the 
directions of these forces are tangential to the shell, the resultant total force is radially 
outward pointing and its size is AF tot = 4AFsin (A</>/2). Correspondingly, AF tot tends 
to 2AFA0 in the A0 — > limit which, in turn, justifies that the force per unit mass, 
due to the non-zero pressure, is limA0->.o AF tot /AM = 87iPR/M she n. Accordingly, the 
equation of motion is of the form 

~ 2M C + M shcll SttPR 

R = + (2 ' 28) 

and the initial condition to this differential equation consists of M s h e ii, Ro and Vo, where 
Vq = dR/dT at Rq. In addition, as in the relativistic case, we need the central mass M c 
as an environmental variable. Formally the gravitational mass can be calculated from 
the initial conditions just as in the Einstein theory, but, as we noted above, in this case 
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the gravitational mass has to coincide simply with M s h e ii- In comparing results relevant 
for the Newtonian and Einstein theories we need to harmonize initial conditions. Note 
that while ro and Ro, or Vq and Vq, have identical interpretations in both theories, in 
setting the value of M she n we have the following inequivalent two choices. We may 
identify M s h e ii with either the rest mass or the gravitational mass in the Einsteinian 
setup. As in the Newtonian regime, i.e. whenever r 3> m c ,m g and Vq <C 1, the 
relation m g ~ m r holds, we could identify either of them with M she n. Nevertheless, as 
the gravitational mass plays the same role conceptionally in both theories, it seems to 
be more appropriate to assume the equality of the two gravitational masses by setting 
M she ii = m g . 

Assuming that P is a given function of the radius P = P{R), the equation of 
motion, ( 12.281) . can be integrated. The existence of such a function is always guaranteed 
whenever an EOS of the form P = P(E) is known because £ = M s h e ii/(47ri? 2 ) itself is 
a function of the radius. It is worth emphasizing that the EOS was kept as completely 
generic throughout the above discussion. 

Multiplying (12.281) by R and by then integrating with respect to T one gets, as the 
analog of (12. lip , the Newtonian energy balance equation 



where W(R) = P(i?)/£(i?). In particular, in the case of a shell with a linear EOS 



where W(R) = W = const, the last term on the right-hand side of (12.291) takes the form 
4W\n(R/R ). 



3. The evolution of multi-layer shell systems 

This section provides a brief review of the analytic and numerical setup we have applied 
in studying the evolution of multi-layer shell systems. This evolution, in general, involves 
a large number of collisions of the shells. 

3.1. Collision of two shells 

As it was emphasized in section IV of [61 j . the description of a collision of two shells 
cannot be given without invoking some further assumptions concerning the interaction of 
the shells. More concretely, there is an ambiguity in the evolution of colliding shells even 
though the energy and momentum conservations are guaranteed to hold as — depending 
on the type of interactions of the shells — more than two shells or even a continuous 
spread of the matter of the original two shells into a thick shell might develop during 
the collision. 

This ambiguity is eliminated if the two shells pass through each other either without 
any interaction, in which case the collision is said to be totally transparent, or when the 
interaction is extremely strong and the ingoing shells merge into a single outgoing shell, 
in which case the collision is referred to as totally inelastic. These two extreme cases are 




(2.29) 



Probing the stability of gravastars by dropping dust shells onto them 

schematically represented in figure El In this paper only the totally transparent shell 
crossings will be considered. The EOS of the shells will also be assumed to be intact in 
collisions. Note that while for dust shells this assumption seems to be appropriate it is 
much less adequate whenever the shells are comprised of strongly interacting particles. 
Below, the basic equations relevant for totally transparent collisions are recalled. 



Shell 4 Shell 3 Shell 3 




Shell 1 Shell 2 Shell 1 Shell 2 



Figure 3. Schematic spacetime diagrams representing the totally transparent (left) 
and totally inelastic (right) collisions. The vertical direction is temporal and time 
progresses upward. 



Detailed investigation of the dynamics of the totally transparent collision of two 
shells can be found in [39]. For this type of collision it is assumed that m r3 = m r i 
and m r4 = m r2 and — provided that the four-velocity of each shell is continuous at the 
location of collision — the momentum conservation can be shown to be equivalent to the 
relations (39] 

p 3 =p x + Ap, (3.1) 
p 4 = p 2 + Ap, (3.2) 

where pi = m r i(dr/dr)i stands for the 3-momenta of the zth shell — here the indexing of 
the shells follows the notation applied on the left-hand side of figure [3] — and 

m r im r2 _ (m gl -h 1 )p 2 -(m g2 + h 2 )p 1 , , 

r c r c - 2(m c i + m gl ) 

where r c is the radius at the collision, m r j and m g j stand for the rest mass and the 
gravitational mass of the ith shell, m c \ denotes the central mass, while w", and 
hi = m r ? I (2r c ) denote the four- velocity, the unit normal and the 'self-gravity' of the ith 
shell, respectively. 

Note that in (13 .ip and (I3.2p both of the positive signs in front of the term Ap are 
correct as, in order to suit to the three-momentum conservation, Ap itself is negative. 
The apparent conflict may be resolved immediately if one takes into account that p 3 

6 A detailed analytic description of the totally inelastic case has been given in the appendix of [9 j while 
the C++ code [ST] is developed so that it is capable of investigating the evolution of shell systems when 
the collisions are assumed to be totally inelastic. 
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and p 4 , as well as p\ and p 2} are components of four- momentum vectors with respect 
to different bases [39] . It also follows from Ap < that, in order to have a transparent 
collision the vector fields ul and have to be arranged so that the contraction u\ri2a 
be positive. 

From the four-momentum conservation the 'energy balance' relations 

m g3 = m gl — AS , (3.4) 

m g4 = m g2 + AS , (3.5) 

can also be derived, where m g i = M 2 — M 1 , m g2 = M 3 — M 2 , m g3 = M 3 — M 4 , 
m g 4 = M4 — Mi; the stands for the mass parameter of the ith Schwarzschild region 
as indicated on the left panel of figure [3] and the measure of the 'energy transfer', AS, 
can be given as 

AS = u\u 2a ■ (3.6) 

r c 

Since both «" and u% are future-directed timelike vectors the contraction ufu 2a is 
negative implying that AS is always positive. This, in particular, means that the outer 
shell loses energy, while the energy of the inner shell increases. 

In subsection 14.31 the above relations will be economized in characterizing the mass 
inflation phenomenon. 

Note that in the case of multi-layer shell systems one should take into account the 
possibility of simultaneous collisions of more than two shells. Since the occurrence of 
these type of events is practically zero, by choosing the time steps to be sufficiently 
small we could guarantee that in each time step merely pairwise collisions occurred. It 
is also important to note that whenever a simultaneous collision of more then two shells 
is allowed to occur the result may always be determined by decomposing the event into 
pairwise collisions and the result is independent of the order of the pairing process. 

Note, finally, that in the Newtonian case the basic equations for totally transparent 
collisions are simple as the masses and the velocities are interchanged in the most obvious 
way. 



3.2. Dynamics of multi-layer shell systems 

In describing the relative motion of multi-layer shell systems we need to specify a 
reference parameter. As long as considerations are restricted to two shells, the most 
convenient reference parameter is the Schwarzschild time coordinate defined in the 
intermediate region while we remain in the region of outer communication, whereas 
inside the black hole region the Eddington-Finkelstein null coordinate serves the same 
purpose. Let us mention that without choosing an appropriate reference parameter one 
may create significant confusion even in the (simplest possible) case of two shells. For 
instance, in jlQ] where the respective proper times measured along the separate shells — 
these, in virtue of (12.23}) and ( 12.24(1 . may differ considerably — were assured to coincide 
and were used as a reference parameter, false conclusions were derived concerning, e.g., 
the extent of regions where the crossing of shells may occur. 
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In describing the motion of multi-layer systems consisting of N (> 2) shells, labeled 
by the index i 6 {1, N}, the Schwarzschild time coordinates defined in the respective 
intermediate regions, associated with any pair of shells next to each other, are — as 
long as only the two shells are concerned — as suitable as before^. However, as was 
discussed before, these Schwarzschild time coordinates cannot be matched properly 
due to their discontinuities across the shells. To overcome this technical difficulty 
the following synchronization process may be applied. In each of the intermediate 
regions the constant t-lines determine simultaneity. The entire spacetime-built up 
from piecewise Schwarzschild regions — is synchronized by applying this synchronization 
process to succeeding intermediate regions starting from the innermost region. 

By making use of this synchronization method, the location, i.e. the coordinates 
of all the shells, can be given as a function of the time coordinate of the innermost 
region, playing the role of reference parameter, t r . In most of the following figures 
instead of plotting the r^\t r ) functions, the expressions Ar^(t r ) = r^\t r ) — f(t r ), 
where f(t r ) denotes the radial center of mass of the r^ l \t r ) distribution, i.e. f(t r ) = 
YliLi Wr l VW(t r )/ YliLi nTv\ will be plotted^. Interestingly, in certain cases, e.g. in 
case of a collapsing shell system where f = f(t r ) is guaranteed to be a monotonically 
decreasing function, instead of t r the radial center of mass f can also be used as a 
reference parameter. For instance, in some of the figures (see, e.g., figures IU [5j [7] and 
E]) the functions Ar^' = Ar^(f) will be plotted. 

It is important to keep in mind that the above-defined synchronization has a certain 
degree of ambiguity concerning the specific i r -value at a given point. Nevertheless, once 
the synchronization is fixed, the crossings of shells are well-defined as the t r -labels are 
uniquely determined, i.e. the motion of the shells may be properly represented in the 
(t r ,r) local coordinates. Note also that using the Schwarzschild time coordinates in the 
intermediate regions we may follow the motion of the shells only up to the appearance of 
the event horizon in either of these regions. One could argue that an external observer 
cannot see what happens to the shells beyond an the event horizon, which manifest 
itself by the infinite value of the Schwarzschild time coordinate. Nevertheless, in some 
cases, such as in the study of mass inflation discussed in subsection I4.3[ it is important 
to descend into the black hole region. Fortunately, just like in the case of a single 
shell, the use of the Eddington-Finkelstein null coordinates instead of the Schwarzschild 
time coordinates resolves the corresponding problem of synchronization for multi-layer 
systems. 

7 If the motion is intended to be described inside the event horizon in these intermediate regions 
Eddington-Finkelstein null coordinates have to be applied instead of the Schwarzschild time as it has 
been indicated several times. Note also that, because of spherical symmetry, we frequently replace the 
spacetime by its factor space with respect to the group SO (3). 

8 By simply reversing the orientation of the synchronization, i.e. by starting from the outside we may 
also use the Schwarzschild time coordinate of the external region as our reference time. 
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3. 3. Initial data for systems 

In specifying the initial data for a shell system, just like in the simple shell case, we need 
as an environmental parameter the central (Schwarzschild) mass, mg, characterizing the 
innermost region. The rest of the initial data set consists of the triples m['', Vq for 
the individual shells, with i G {1, 2, N}, which are required to satisfy the relations 
?o < To for i < j. In determining the gravitational mass, rrig , of the ith shell — by 
making use of the pertinent form of (12.16!) — we nee d to know the central mass, m['' , 
felt by the ith shell. It can be justified that for % > 2 the relation mc = mi* -1 ^ + m[ ! ^ 
holds, with mc = wig. Note that the initial data m.Q , \ also have to satisfy the 
pertinent form of (12. 17j) for each individual shell. 

Note also that, in specifying the initial state of our multi-layer shell system the 
EOS of the individual shells have to be fixed. 

3.4- The numeric algorithm 

The open source C++ code of the applied numerical algorithm can be downloaded (see 
|57J). This subsection is to provide a short outline of this numerical algorithm. 

Apart from shell collisions, the motion of the individual shells can be treated 
separately when we suppress the indices of shells hereafter. In order to be able to 
integrate the equations of motion — see (12. lip — we need to know m c , m g and m r (r). 
Once m c and suitable initial data — consisting of mo, r and t>o — are specified the value 
of m g and m r (r) can determine as on one hand m g is given by (12. 16ft while on the other 
hand, m r (r) = 4ira(r)r 2 and the functional relation a = o~(r) can be deduced with the 
help of an EOS, as was described at the end of subsection 12.11 For instance, in the 
particular case of a dust or fluid shell with a homogeneous linear EOS, m r (r) can be 
given analytically by making use of (12.150 . In all the other more generic cases numerical 
algorithms can always be used in solving (I2.14p . 

Once m r (r) is known, the absolute value of dr/dr, along with (dt±/dr), can be 
determined as a function of the radius by making use of (12.111) . (12.231) and (12.241) . The 
sign in front of the velocity dr/dr depends on the initial data and its value remains 
fixed until the appearance of a turning point of the shell if it occurs at all. The velocity, 
with respect to the Schwarzschild times dr/dt±, can also be given by applying the chain 
rule 



In describing the motion of a multilayer shell system, proceeding from the inside to the 
outside, we need to apply the synchronization outlined above. Therefore, the motion 
of each shell is determined with respect to both the inner and outer Schwarzschild 
time (or the Eddington-Finkelstein null coordinates) because — as follows from the 
synchronization process — in specifying the 'time step' in the outer regions we need to 
know the 'time step' applied in the innermost region. 




(3.7) 
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In integrating the equations of motion the fourth-order Runge-Kutta algorithm was 
applied. The appearance of a turning point in the integration process is indicated by 
the fact that the right-hand side of (12.111) becames negative. When this happens the 
code does not undertake the last step and by making use of a root finding the location 
of the turning point is determined. The equation of motion (12. lip is then integrated 
such that the sign of the velocity, dr/dr, calculated from it is changed to the opposite 
of its previous value. A change in the order of the shells, i.e. the ordering of their radial 
coordinates, always indicates that a collision occurred somewhere. Then the evolution of 
the system is then held on, for a short while, until the precise location (in space and time) 
of the collision is determined by linear approximation. At that event, by making use 
of the conservations equations and the assumption about the EOS, we determined the 
initial data for the new shell(s). (If inelastic collisions are allowed to occur the number 
of shells has to be decreased by 1.) The evolution of the full multi-layer system is carried 
on with the inclusion of the yielded new shell(s) afterwards. Once the innermost shell 
reaches the origin, the evolution of the system is continued so that it is taken out of the 
system of shells and its gravitational mass is added to the center mass of the innermost 
region. Whenever the Eddington-Finkelstein coordinates were applied, the algorithm 
was basically the same with the distinction that derivatives dt±/dr were replaced with 
the corresponding expressions dv ± /dr. 

The precision of the applied numerical schema is checked in the case of a system 
formed by two repeatedly intersecting equal rest mass dust shells (see figure [11] in 
section |4~3]) . This justifies that our numerical code is convergent even though succeeding 
collisions occur. 

4. The main results 

One of the most interesting questions in working with systems of thin shells is how to 
mimick thick shells. However tempting certain analogies might be, one should keep 
in mind that even though one is using a huge number of thin shells, the continuum 
distributions cannot be properly modeled, with the exception of the dust case, because in 
the thin shell limit the interaction of particles in the direction transversal to the shells — 
apart from the gravitational one — is neglected. What is possible to do consistently is to 
investigate an approximate model using a large number of dust shells. Let us mention 
here that the only attempt to investigate such an approximate model — as far as we 
know — was made in [55] although the number of the shells was kept minimal, i.e. only 
the evolution of a two-shell system was investigated. In contrast to this, with the help 
of our C++ code [57], the evolution of a large number of shells can be determined. 

4-1. The study of 'simple 7 systems 

In order to provide a better understanding of the dynamics of our thick shell mimicking 
the multi-layer thin shell system as reference solutions, the evolution of N = 16 shells 
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with uniform initial mass and radial distributions will first be examined. Since even 
in this case there are too many parameters characterizing the system, it seems to be 
appropriate — as will be done below — to proceed in small steps by changing the initial 
conditions of the reference solution almost parameter by parameter. 

As mentioned above we shall consider the time evolution of shell systems consisting 
of iV = 16 shells, with the exception of panel (c) in figure H] where the evolution of 
N = 64 shells will be considered. The total rest mass YliLi m o^ °f these shells will be 
72, while m§ = 0, in each of the following cases. The mass distribution will be uniform 
on each panel of figure HI while it is uniform on panel (e) and centered on panel (f), 
possessing the mass distribution as specified by (14. ip . The initial radial distribution j-g , 
in the case of uniform distributions, is chosen so that ?q = 1999 + i, where i runs from 1 
to 16. It is centered on panel (f) and it is random on panels (g) and (h) in figured! It is 
important to keep in mind that even in the case of a uniform mass distribution, only the 
rest masses m!^ could be arranged to be equal to each other, whereas the corresponding 
surface mass densities differ slightly according to the relations m = 47T(7qVo^. 
For simplicity, only dust shells are considered — with the exception of panel (d) in figure 
0] — and the initial velocity, Vq , was chosen to be zero, i.e. all the shells start from rest. 

On each of the panels of figures H] and [5] on the horizontal axis the radial center of 
mass f is indicated, which means that the synchronized events are arranged so that they 
lie along vertical segments. The initial radius distribution of the shells is indicated by 
strokes on the right side of each plot. In order to help the recognition of the developing 
structures the following type of gray shadowing had been applied. A lighter gray 
shade was used between the momentary outermost and the last but one shell, while 
the complementary inner part received a darker gray coloring. Since we intended to 
compare the Newtonian and fully relativistic time evolutions in the relativistic case, all 
the plots were produced by making use of the Schwarzschild time coordinates, whence 
all the evolution stops before reaching 'the even horizon'. 

Let us provide now a brief outline of the plots in figures 0] and [5j The letters applied 
in the following itemization refer to the labels of the the individual panels. 

(a) In the Newtonian setup the time evolution of the above-specified shell system with 
an uniform radius and mass distribution is shown. 

(b) The time evolution, with initial data corresponding to that of panel (a), is depicted 
in the fully relativistic case. 

(c) This justifies that the basic characters do not change drastically whenever the 
numbers of the shells is increased by a factor of 4. The time evolution of 64 shells 
is depicted in the fully relativistic case. 

(d) The fully relativistic time evolution of 16 shells is considered using the same initial 
configuration as in panel (b) with the distinction that the EOS of the shells is 
homogeneous linear, V = wa with w = 0.006.§ 

9 The order of magnitude of the value of w was determined as if the shell system formed a virialized 
system in the kinetic gas theory. 
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(a) uniform distribution of 16 shells (newtonian dynamics) 
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(b) uniform distribution of 16 shells (relativistic dynamics) 
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(c) uniform distribution of 64 shells (relativists dynamics) 
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(d) linear equation of state with w=0.0025 (relativistic dynamics) 
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Figure 4. The time evolution of shell systems with uniform mass and initial radial 
distributions is shown. On the horizontal axis the radial center of mass, f, of the radius 
of the shell system is indicated, while on the vertical axis the deviation Ar^ = r^ % ' — r 
relevant for the individual shells is plotted. The upper two panels are to indicate the 
similarities and differences between the Newtonian and fully relativistic evolutions. 
Although it might look strange at first, it is important to keep in mind that for the 
depicted collapsing systems — where f is monotonically decreasing — here, and in figures 
EJ Ul HI M andflOl time progresses from right to left. 
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(e) centralized mass distribution of 16 shells (relativistic dynamics) 
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Figure 5. The time evolution of shell systems are shown. In the top panel the mass 
distribution is uniform while in the other three panels the initial radial distribution 
is non-uniform; in fact, the initial radial distribution is chosen to be random on the 
lower two panels. As in figure [4] the radial center of mass, r, of the radius of the shell 
system is indicated on the horizontal axis, while on the vertical axis the deviation 
Ar^ = r-W — f relevant for the individual shells is plotted. In spite of the significant 
differences in the initial part of the indicated evolutions, the similarities in the final 
parts are considerable. 
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(e) In this plot the evolution of a system with a uniform initial radial distribution but 
with the centralized 

m if 1 < z < 8; , , , . 

m 1 ' = < ' , ~ (4.1) 
1 16-i, if 9 < % < 16. v ' 

mass distribution is shown. The latter is indicated by the thickness of the horizontal 
strokes at the right edge of each of the individual figures. 

(f) In this panel the evolution of a system, with uniform initial mass distribution 
and with symmetrically centered radial distribution, is shown, where 7g = 2000, 
r{, 16) = 2015 and 

rf +1) - rf = (9 - i) ■ A , for i = 1, 2, 7 ; (4.2) 
r o ] - r o' 1] = (* - 8 ) • A > for i = 10, 11,... ,16, (4.3) 
with A = 15/71. 

),(h) In these two panels the mass distribution is uniform while the initial radial 
distributions are uniform random distributions, although the initial location of the 
outermost shells were chosen to be the same as on all the other plots, i.e. = 2000 
and rj 16) = 2015. 

By comparing the plots in figures H] and [5] the following conclusions may be drawn. It 
is clearly visible that whenever the shells are starting from rest, not only the nearby shells 
cross but the crossing of all the shells seems to be generic. It has been justified by our 
numerical experiences-see also panel (c) in figure @] — that the developing basic structures 
seem to be independent of the number of shells. Let us mention that shell crossings are 
known to occur even in the continuous model of spherical dust systems [63J. Note, 
however, that they are also frequently referred to as 'shell crossing singularities'. These 
type of singularities are known to be much weaker than the central ones; nevertheless, 
in the continuous model the evolution stops whenever they occur. 

Returning to the interpretation of panels (a)-(d) of figure HI and (g)-(h) of figure El 
by comparing these plots it is straightforward to recognize that either type of deviation 
from the uniform distribution yields a visible dispersion of the shells. The focusing of 
the shells is lost faster and the spreading is more intensive as the initial distributions 
get further and further away from uniformity. 

Let us now make some more specific comments. Comparing panels (a) and (b), 
depicting the Newtonian and relativistic evolution of completely uniform distributions, 
it is visible that the nonlinearity of the relativistic evolution yields a larger dispersion 
indicated by the fact that in the Newtonian case the first five knots are pretty well 
focused, while the loss of the coherence starts earlier in the relativistic case. Panels 
(c) and (d) of figure EB justify that the main characters of the evolution do not change 
whenever either the number of the shells is increased or the EOS is changed a little. 
In panel (e), where the initial mass distribution is centered, as it is specified by (14. ip , 
knots do not develop and at the beginning the dispersion of the system dominates. 
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Figure 6. The time dependence of the differences of the radial center of mass values 
of the shell systems depicted in panels (a), (b) and (e) in figures @] and O and that of 
the Newtonian reference solution of a single shell with mass M s h e ii = 16 are indicated 
by plotting the f(t c ) — r-^s(t c ) functions, where t c denotes the Minkowski time in the 
central region. f(t c = 0) = 2007.5 and f — r^s = at t c = for each of the individual 
systems. 



On the other hand, during the second half of the indicated period the shells start to 
be concentrated around two of the highest mass shells, and then the groups formed 
this way start to oscillate around each other. In panel (f), where instead of the mass 
distribution the initial radial distribution is concentrated, a similar pair of groups of 
shells can be seen to develop. In both panels (e) and (f) it is also visible that the radius 
of the outermost shell varies on a significantly larger scale than in the other panels. 
In particular, their oscillation amplitudes increase during the initial part and start to 
decrease only later. Nevertheless, the relative variety on their motion, with respect to 
the characteristic size of the shells composed of all but the outermost shells, remains 
significant during the entire evolution. The gray shading makes this type of effect more 
apparent. In panels (g) and (h) the initial part of the distribution of the shells is random. 
It is clearly visible in these panels that the shells are not concentrated into two groups 
as before. 

In order to make the differences of the above-discussed dynamics more 
approachable, the time dependence of the differences, f — tns of the radial center of 
mass of the shell system depicted in panels (a), (b) and (e) in figures H] and |5] is shown 
in figure O where t c denotes the Minkowski time of the central region and r^s = ^ns(^c) 
stands for the time dependence of the Newtonian reference solution of a single shell 
with mass M s h e ii = 16. In the left panel of figure |6] a zoom into the final part of the 
evolution is shown in order to make the slight differences more visible. Note that the 
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wavy figures do correspond to the formation knots indicated by the initial parts of the 
time evolution in panels (b) and (e) of figures H] and [5j It is worth emphasizing that the 
radial center of mass, f = f(t c ), is a monotonically decreasing function for each of the 
individual systems. 

4-2. Dispersion of small perturbations 

In this subsection the number of shells is increased significantly from N=16 to 101 in 
our reference simulation with uniform initial mass and radial distributions shown by the 
upper-left panel of figure [7J 

The applied small perturbations on the evolution of equal mass shell systems are 
produced either by doubling the mass of a single shell (or ten shells) of an initially 
uniform distribution, or by taking out a single shell (or ten shells) from an initially 
uniform distribution. The initial part of the evolutions of these slightly perturbed 
systems are compared to that of an initially uniform distribution in figure [7J 

The evolutions relevant for systems where the rest mass of one or ten shells is 
doubled are shown in the middle-left and lower-left panels respectively. The change of 
the motion of the other shells is noticeable but at first glance it is not too striking. 
As opposed to this, in the panels on the right column of figure UJ — where in the upper 
and middle panels only a single shell is removed, while ten shells are taken out from the 
initially uniform distribution in the lower panel — the perturbations yield more significant 
changes. The consequences of the indicated small changes are more significant in spite 
of the fact that the change of the total rest mass is about 1% whenever only a single 
shell is left out from the reference simulation with uniform initial distributions. 

By the inspection of the figure above, the following important qualitative 
observations can be made. First of all, there is a tendency for the formation of a 'crust'— 
represented by the increase of the density of shells — at the edges of the widening gaps. 
Second, following the widening of a gap, a reversing of the sides also occurs, i.e., the 
innermost shell becomes the outermost and vice versa, as it is clearly visible in the 
colored figure [H In addition, the crusts are 'growing' and during the collapse the entire 
system gradually becomes a dispersed two-shell system. Third, around the anti-gaps an 
increase in the density of shells is also noticeable, although the growing rate is much 
lower than in case of gaps. The two lower panels in figure [7J indicate that the evolution 
of systems with gaps and anti-gaps get closer to each other when a large number of gaps 
and anti-gaps are uniformly distributed in the initial configurations. 

4-3. Mass inflation 

In starting this subsection we would like to mention that some of the arguments below 
were motivated by claims of [39] about the time evolution of a system formed by a pair of 
repeatedly intersecting equal rest mass shells. However, we would like to emphasize that 
our conclusions about the possible rate of the blowing up of the mass of the intermediate 
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Figure 7. In these plots the initial parts of the evolution of shell systems yielded 
by slight perturbations of a reference system — consisting of 101 uniformly distributed 
equal mass shells starting from rest, depicted by the upper-left panel — are shown. The 
initial radial distribution, r^ 1 , of the reference system is chosen such that the location 
of the shells is symmetric to with r$ = 9999 + i, while mj' = 1 and = for 
i = 1,2, ...,101. The systems with 'gaps' shown on the upper-right and middle-right 
panels are the result of removing the 51th and 31th shells from the reference system, 
respectively. The middle-left panel shows a system with an 'anti-gap' where the mass 
of the 31th shell is doubled, i.e. mj 31 ' = 2. In the lower panels the evolution of systems 
with ten uniformly distributed gaps or anti-gaps is shown. The radial center of mass f 
is indicated on the horizontal axis — which takes values from the interval [9425, 10050] 
on each panel — while in the vertical direction the deviations, Ar'*' = — f, are 
shown. 
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Figure 8. In this colored plot a longer period of the time evolution, for the system 
in the upper-right panel of figure [7j is shown. It demonstrates that a slight change of 
the uniform initial data might result significant differences in the long run, and the 
formation of 'crusts' and reversing of sides are also noticeable. 



regions — these are also justified by means of numerical investigations below — differs from 
that of [39]. 

In describing the unbounded growth of mass, recall first that the squares of energy 
and momentum exchanges, AS and Ap — for their definitions see subsection 13 .11 — are 
related as 

AS 2 - Ap 2 = (4.4) 



r, 



c 



which implies that for AS the inequality 

AS > (4.5) 

r c 

holds. Now, by making use of the relations m g 2 = M3 — M 2 and m g 4 = M4 — Mi, along 
with the vanishing of m c \ = Mi, in virtue of (13.51) . we get 

M 4 + M 2 = M 3 + AS . (4.6) 

Hence, in virtue of H4.5[) and (14.61) . for the radial center of mass, M = (M4 + M 2 )/2, of 
the mass parameters of the successive intermediate regions, the relation 

M>Um 3 + ^^) (4.7) 



2 V r c 

has to hold. This latter inequality, whenever r c tends to zero and AS becomes much 
larger than M 3 — for the case of colliding shells possessing linear EOS with wi = and 
W2 = c 2 2 — implies that M must tend to infinity such that the asymptotic blow-up rule 

M > Cr- (1+2wi+2w2) (4.8) 

holds. 

The process, which in repeated collisions of a pair of shells leads to an unbounded 
growth of the mass parameter of the intermediate region, is the mass inflation. The blow- 
up behavior of some simple configurations, along with the justification of the estimate 
(14. 8 p , is shown in figure M 
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crossing radius 

Figure 9. The log-log plot of the mass parameter of the intermediate region with 
respect to the radial coordinate of the collisions is shown for a system formed by two 
repeatedly intersecting equal rest mass shells possessing linear EOS with w% and W2- 
Note that the dots in this figure signify, at the location of the collision, the value of 
the mass of the intermediate region produced in the collision, while between any two 
collisions the mass of the intermediate region is constant. The plot with w\ = = 
corresponds to the case of colliding dust shells. In both cases the initial data were 
synchronized in the Schwarzschild time of the intermediate region and was chosen to 
be the same as in [55], i.e. such that mg — 0, mj 1 ' = 100, mg' 1 ' = 100, mj 2 ' = 100, 
m g ( 2 ) = 90.05906. (Note that in virtue of (|2.16j) m s ^ and m g ( 2 ) determine the initial 
values Wq 1 ^ and v^\) The low radius behavior of the plotted curves justifies that the 
Schwarzschild mass of the intermediate regions grows as M ~ r~ a , where for a the 
fits to the estimate a = 1 + 2w\ + 2u> 2 formulated by (|4.8|) . The arrows point to the 
locations where the critical radii of the external region become larger than the radii of 
the succeeding collisions, i.e. where e t+ for the outer shells change their signs. 

Some remarks are now in order. In explaining the relatively small values of w it 
should be mentioned that whenever larger values of u>i and w 2 are applied, or W\ and 
W2 differ significantly, the fluid shells start to move outwards and it may happen that 
they collide only once or they do not collide at all. It is also important to note that 
whenever the collapse and mass inflation occur the blow-up rate behavior is insensitive 
to the initial data of the shells. 

Note that mass inflation is not new; it is known to occur in the continuum limit (see, 
e.g, |64j for a recent numerical investigation). Nevertheless, at first sight the occurrence 
of the mass inflation in the thin shell formalism is unexpected because it is known that 
the mass of the intermediate region cannot be larger than that of the outer region unless 
the radius of the collision, r c , becomes smaller than the critical radius, r, of the outer 
shell. Note that such a critical value in the case of the dynamics of a single shell is 
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Figure 10. The values of finner and ^outer as a function of r c are plotted on a log-log 
scale for the case of repeatedly colliding equal rest mass dust shells corresponding to the 
'points' represented by the '+' sign in figure [9l The intersection of the r c = ro U ter(?"c) 
curves signifies — indicated by the arrow — the location where e t+ of the outer shell 
changes sign. Note also that gravitational mass of the outer shell changes sign where 
^outer attains its maximum. 



extremely small, much smaller than the Schwarzschild radius. In order to show that the 
above-formulated expectation is justified by the investigated time evolutions in figure 
ITUl the time dependence of the critical radii of the inner and outer shells — their relative 
location varies in time — along with the time dependence of the radius of the collision 
is plotted. The location where for the first time r c becomes smaller than r for the 
temporarily outer shell is the very location where the mass of the intermediate region, 
M4, becomes larger than the mass of the outer spacetime, M 3 . 

The precision of the applied numerical schema is determined in the case of the 
above-described system formed by two repeatedly intersecting equal rest mass dust 
shells. Denote by the numerical value of the radius of the (temporarily) internal shell 
relevant for resolution A. In figure fTTl the time dependence of the difference r^ ^ — r^j. 
is plotted for the initial period for A = 10 n 5 with n = 1,2,3 and r-^j denotes the 
reference numerical solution with the smallest resolution 5. As is expected for any fixed 
value of n, the difference rj£° ^ — rj^ is increasing in time, however, by decreasing the 
value of n by 1 yields an order of magnitude downward shift of the successive curves. 
This justifies that our numerical code is convergent even though collisions — indicated 
by the jumps on the curves — occur. 

Let us finally mention that mass inflation is not specific only to two-shell systems as 
it did occur for systems consisting of more than two shells. According to our experience 
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Figure 11. The time dependence of the difference, S>> — r^, of the numerical 
value of the radius of the (temporarily) internal shell is plotted relevant to resolutions 
10 n S with n = 1, 2, 3 and for the reference numerical solution r^. 

the masses of each of the intermediate regions increase apparently following the above- 
derived power law rule, but there was a definite order kept during the evolution. 
More definitely, at any instant of the Eddington-Finkelstein time the mass of an inner 
intermediate region was, in all of our simulations, smaller than the mass of any of the 
intermediate regions located outwards with respect to the inner one. 

5. Final remarks 

The relativistic time evolution of multi-layer spherically symmetric shell systems has 
been investigated. After recalling the basics of the analytic setup a newly developed 
numerical code is introduced. This numerical method was made to be capable of 
following the time evolution of systems comprising of great numbers of colliding shells 
such that whenever collisions occur they are assumed to be totally transparent. 

By making use of our numerical method for the first time, the relativistic time 
evolution of numerous shell systems involving large number of thin shells could be made. 
As these systems may be considered as approximate models of thick shells, the results 
reported in this paper provide insights about their dynamical behavior as well. The 
most important observations we have made can be characterized by the key phrases: 
concentrations of subsets of shells, formation of 'crusts' at the boundaries, reversing of 
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the sides of gaps. 

The chosen analytic setup ensured that the evolution of the considered shell systems 
can be investigated both in the domain of outer communication and in the black hole 
region. This made our numerical code capable of studying mass inflation within the thin 
shell formalism. We would like to emphasize that beside the numerical investigation of 
mass inflation, an estimate explaining the main features of the blow-up behavior of the 
mass parameter of the intermediate region is also provided. 

As was mentioned in the introduction, there are a great number of astrophysical 
systems which can be modeled by shells. For instance, there are plans to investigate 
the interaction of repeated quasi-spherical matter ejections by supernovas using the 
developed method. 
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